% fig1a_reproduction_fsolve.m — 用 fsolve 改进后的 Fig.1(a) 复刻
clear; clc; close all;

%% —— 参数设置 —— 
g11 = 1; g12 = 1; g21 = -2; g22 = 2;    % 连接权重
mu_vec = linspace(0,1,200);             % μ 从 0 到 1，共 200 个采样点
nMu    = numel(mu_vec);

% fsolve 选项
opts = optimoptions('fsolve', ...
    'Display','off', ...
    'FunctionTolerance',1e-8, ...
    'OptimalityTolerance',1e-8, ...
    'StepTolerance',1e-8);

%% —— 同时求解 (x,y) 定点 —— 
X_all = cell(nMu,1);  % 存储每个 μ 的所有 x 定点
for i = 1:nMu
    mu = mu_vec(i);
    sols = [];
    % 在 [-3,3]×[-3,3] 区域，用 7×7 网格做初值
    grid1 = linspace(-3,3,7);
    for x0 = grid1
        for y0 = grid1
            try
                sol = fsolve(@(uv)[...
                    mu*uv(1) + g11*sin(uv(1)) + g12*sin(uv(2)) - uv(1); ...
                    mu*uv(2) + g21*sin(uv(1)) + g22*sin(uv(2)) - uv(2)], ...
                  [x0; y0], opts);
                % 只保留实部，且抛弃明显发散的解
                sol = real(sol);
                if all(abs(sol)<20)  % 合理范围之内
                    sols(end+1,:) = sol.'; %#ok<AGROW>
                end
            catch
                % 无收敛则忽略
            end
        end
    end
    % 去重（按 x,y 同时去重）
    if ~isempty(sols)
        sols = unique(round(sols,4),'rows');
        X_all{i} = sols(:,1);  % 只取 x 分量
    else
        X_all{i} = [];
    end
end

%% —— 绘制定点分布 —— 
figure('Position',[100,100,700,500]); hold on;
for i = 1:nMu
    mu = mu_vec(i);
    Xs = X_all{i};
    plot(mu*ones(size(Xs)), Xs, '.k', 'MarkerSize', 6);
end

%% —— 四组初始值 & 颜色 —— 
init_vals = [ 0.1,  0.1;
             -0.1, -0.1;
              0.1, -0.1;
             -0.1,  0.1 ];
colors = [1,0,1;      % magenta
          0.6,0.3,0;  % brown
          0.5,1,0.5;  % prasinous
          0,1,1];     % cyan
labels = {'(0.1,0.1)','(-0.1,-0.1)','(0.1,-0.1)','(-0.1,0.1)'};

%% —— 迭代演化并绘制轨迹 —— 
Niter  = 500;   % 总迭代步数
Ntrans = 400;   % 舍弃前 400 步
for i = 1:nMu
    mu = mu_vec(i);
    for j = 1:4
        x = init_vals(j,1);
        y = init_vals(j,2);
        for k = 1:Niter
            x_next = mu*x + g11*sin(x) + g12*sin(y);
            y_next = mu*y + g21*sin(x) + g22*sin(y);
            x = x_next;  y = y_next;
            if k > Ntrans
                plot(mu, x, '.', 'Color', colors(j,:), 'MarkerSize', 4);
            end
        end
    end
end

%% —— 美化 —— 
xlabel('\mu','FontSize',12);
ylabel('x','FontSize',12);
title('DBHNN 模型 \mu–x 平面定点与并存分岔轨迹','FontSize
